{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set some ipython notebook properties\n",
    "%matplotlib inline\n",
    "\n",
    "# set degree of verbosity (adapt to INFO for more verbose output)\n",
    "import logging\n",
    "logging.basicConfig(level=logging.WARNING)\n",
    "\n",
    "# set figure sizes\n",
    "import pylab\n",
    "pylab.rcParams['figure.figsize'] = (10.0, 8.0)\n",
    "\n",
    "# set display width for pandas data frames\n",
    "import pandas as pd\n",
    "pd.set_option('display.width', 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-4.12839474 -0.71069159  7.83908632]\n",
      "[[ 0.17830987  0.86914323  0.46129777]\n",
      " [ 0.65036916 -0.45590689  0.60759268]\n",
      " [-0.7383939  -0.19167408  0.64655665]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "K = np.array([[1., 2, 3],[2,1,5],[3,5,1]])\n",
    "S,U = np.linalg.eigh(K)\n",
    "print S\n",
    "print U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-4.12839474, -0.        ,  0.        ],\n",
       "       [-0.        , -0.71069159,  0.        ],\n",
       "       [-0.        , -0.        ,  7.83908632]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.eye(3) * S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.17830987,  0.86914323,  0.46129777],\n",
       "       [ 0.65036916, -0.45590689,  0.60759268],\n",
       "       [-0.7383939 , -0.19167408,  0.64655665]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 2., 3.],\n",
       "       [2., 1., 5.],\n",
       "       [3., 5., 1.]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(np.dot(U,np.eye(3) * S),U.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.04347826  0.56521739  0.30434783]\n",
      " [ 0.56521739 -0.34782609  0.04347826]\n",
      " [ 0.30434783  0.04347826 -0.13043478]]\n",
      "[[-1.04347826  0.56521739  0.30434783]\n",
      " [ 0.56521739 -0.34782609  0.04347826]\n",
      " [ 0.30434783  0.04347826 -0.13043478]]\n"
     ]
    }
   ],
   "source": [
    "print np.linalg.inv(K)\n",
    "print np.dot(np.dot(U,np.eye(3) / S),U.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.  0.4 0.6]\n",
      " [0.4 1.  1. ]\n",
      " [0.6 1.  1. ]]\n",
      "[[-2.55351296e-15 -5.00000000e+00  5.00000000e+00]\n",
      " [-5.00000000e+00 -1.60000000e+01  1.90000000e+01]\n",
      " [ 5.00000000e+00  1.90000000e+01 -2.10000000e+01]]\n",
      "[[0.59130435 0.11304348 0.06086957]\n",
      " [0.11304348 0.73043478 0.00869565]\n",
      " [0.06086957 0.00869565 0.77391304]]\n"
     ]
    }
   ],
   "source": [
    "V = .2 * K + .8 * np.eye(3)\n",
    "print V\n",
    "print np.linalg.inv(V)\n",
    "print .2 * np.dot(np.dot(U,np.eye(3) / S),U.T) + .8 * np.eye(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.08026756  0.04347826  0.02341137]\n",
      " [ 0.04347826 -0.02675585  0.00334448]\n",
      " [ 0.02341137  0.00334448 -0.01003344]]\n",
      "[[-0.08026756  0.04347826  0.02341137]\n",
      " [ 0.04347826 -0.02675585  0.00334448]\n",
      " [ 0.02341137  0.00334448 -0.01003344]]\n"
     ]
    }
   ],
   "source": [
    "print np.linalg.inv(13 * K)\n",
    "print 1./13 * np.linalg.inv(K)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "print K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "factor = np.diag(K).sum()\n",
    "print factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.33333333 0.66666667 1.        ]\n",
      " [0.66666667 0.33333333 1.66666667]\n",
      " [1.         1.66666667 0.33333333]]\n"
     ]
    }
   ],
   "source": [
    "factor = 3\n",
    "Kn = K/factor\n",
    "print Kn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.49965694678637196"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.std(Kn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1. 1.]\n",
      "[[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "s,u = np.linalg.eigh(np.eye(3))\n",
    "print s\n",
    "print u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 0., 1.]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(np.dot(u,np.eye(3) * s),u.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "u2 = 1\n",
    "print np.dot(u2,np.eye(3)*s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "588.7132803438536 283.0759192455774\n"
     ]
    }
   ],
   "source": [
    "#Why can we assume that h2 and e2 add to one?\n",
    "#Because y is unit standardized then its s.d. is 1\n",
    "# and K's diagional is N\n",
    "\n",
    "covar = np.array([[float(num)] for num in xrange(490)])\n",
    "np.random.seed(0)\n",
    "pheno = covar * 2.0 + 100 + np.random.normal(size=len(covar))*10\n",
    "print np.mean(pheno), np.std(pheno)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAHVCAYAAABWhEeLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xt0lPed5/nPr6okIUBIgDDIIKEQI5uLDQbaydjBYCDBTjxx9mw3znSPx1bSx5uT7pnu9XjTcbY7O864HcPpE7fb3ScJOLODc7KRe3JZX7aTHmMwwjixI4htGYubEyyDKDCXEmAkMKrf/vE8JZVKdXnqkVTX9+scG1QUqKTqJl/X96n3z1hrBQAAgNwJ5PsBAAAAlBsGMAAAgBxjAAMAAMgxBjAAAIAcYwADAADIMQYwAACAHGMAAwAAyDEGMAAAgBxjAAMAAMixUL4fQDr19fW2ubk53w8DAAAgoz179pyy1s7wct+CHsCam5vV0dGR74cBAACQkTHmPa/3ZQUJAACQYwxgAAAAOcYABgAAkGMMYAAAADnGAAYAAJBjDGAAAAA5xgAGAACQYwxgAAAAOcYABgAAkGMMYAAAADnGAAYAAJBjDGAAAAA5xgAGAACQYwxgAAAAOcYABgAAkGMMYAAAADkWyvcDAAAAGC89B7v00+//nXovn9eMj31c/9sD3873Q5LEK2AAAKBE7frxz/X//M3/oUtHT2jCyYs6/9rb+v53Hsr3w5LEK2AAAKAEbf2H7+iD3dslSYFggwIVcxT96KiOHTmc50fmYAADAAAl5dG/+VNVHgxLkoKVi3VVzTrVVwR16qMrOmV25PnRORjAAABASejc0aEXnnlUlWcvS3KGr2vqPqMlE0MykqITgjpbuSS/D9LFAAYAAIpebOVYIWflOGPyp7R4cpPqK4OSJCOjgKwWT+rP7wN1MYABAICilrhyvKbu01oysULGOL9uZGStlTFRTR7Yl8dHOoQBDAAAFKXOHR16vu0xVUacV7WGrRyNM3hJUthEdDx4VhcvTta0qu58PuRBDGAAAKDo/OK7P9a+l3+kSjkrx+CE5Vpcc61aqodWjmETUWfwPb0+dUDH6up1deSYGk/9gT6d34cuiQEMAAAUmScfeVCXOvdLGnqX4w0TQ6oLOa94GRl1BY5qd8UBvdMwV6/Mv0HWGAVtVLe8/XY+H/ogBjAAAFDwYkX7yIWIqiOXJI1cOUpDw9dPp5/Ubxv/QO/VN0gykjEakFVkdmE06BnAAABAQdv145/rtf/3B5Kkag2tHOdNatGSiSEFjCQZnTC9eit4RL+YM0m75n9K1rjDljGStQrI6hr9RtJ/yNNXMoQBDAAAFKzEon2waoHqq2/Qwgkh1VfGXs0avnLc1bJE1n3VS5JkrYyium9gi6oOH8jPF5KAAQwAABSkZEX7xqqgmiuDnlaOslaSFNCAWge26K62HTrSOy1PX81wDGAAAKCgZCraG3flKEldgaPa0tSfdOUoRbXC/kZ3Rp/V8rZ3VfnroLrWVufjSxqBAQwAABQML0V7a61OBJzExC/mTEq5cvzSwGbd0fmSJr4Y1FuhkJ7/44BWXfhY/r64OAxgAACgIHgt2u8PHkt+vVfCyvELbTtU/WqFtqwPaPvSgG7vXqHrohPz9eUNwwAGAADyKpui/aHQce2c2pfieq/hK8cJrwYHh69rTqzU7P56dVw4p8/k70sdxAAGAADyJrFon2zlGF+03z+zUQcamhQ1zq/Hhq9kK8fn/n1Ah2Yb1YTX6beRtXpLUd08+XD+vtg4DGAAACAvEov2yVaOiUX7qDFSiuu9EleOxhqFwrcrGr5Wy/t/q2MTZqnvqr78fLEJGMAAAEBO9Rzs0tOPf1OhM84wlGrlmDYvkeR6r/iVY22kUT29n9N1Jz7U6tPPyiiqARPUh7O5BgwAAJSZWNU+pNRF+9jKMXXRfkDLbIdqFdHK6E4ta3tX+04OrRynhNfqSnix1ve+oXkXj8h9zUwBe0VTT1Xl7WuPxwAGAAByIr5qHwur+inatw5s0R1vvaTgeaOq1wPa2lKh7SuHVo613fVaffpZBWTdP1WSrGSMWuobc/1lJ8UABgAAxlXPwS798ImHFTx1QVL2h2gnWzdO3F2haEB66jMjV463nd4l44xuLivJqH56re784r/N8VefHAMYAAAYN7GVo/uexWHDl9dDtBPzElW/CurFpUbt1w+tHI9G1mnhuXeGDV/O2CVZGU2/ulH3PtAqNd6Uh+/CSAxgAABgXCQepJ2YmPC6chyRl/gTZ/Ay1qgmvFYD4UW6o/eXw673csYuo48m12nRogX6tw/8n7n/BqTBAAYAAMacl6p9yuErQ9F+wcnrVHdqqi5+2KK6E9Hk13vJqH9Wo9bfukAr/tf/lPOvPxMGMAAAMGbSVe1jK0cpQ2IiTdF+bniVOiLrFZUZsXKMsTKqnNWoDfd8QY0r1uf4O+ANAxgAABgTmar2XhITqYr2h68OqCa8Vm9H1mpWf1g3JiQmYq96WUkzr7le9/ztt3P/DcgCAxgAABi1bKv26a73Slw51p+frYrudYqenapVF9q18HyXggkrRyujKwV6vVcyDGAAAMCXnoNden9fp15q/7kqe85L8lm1V+qi/ZTwWv3efYdjrGg/9KqXw8qoasYM3X3XJ9T46a/k9pvgEwMYAADI2q4f/1yvP/sDWavBleNg1X5SyB2SMq8c46/3Sla0HwgvTvoOR+ffzkcz58zVPX95b8EkJrxgAAMAAFnZvHGTzu1tl+Qe8TOsam/c27NMTGwLei7ax9aNNhTS4mvnF8XKMREDGAAA8MQ5RPthhc44RftAsEHBqgWaN3lJ0uu9vK4c46/3Sl+0d171qq2fpuuWXq3mG28r2Hc5ZsIABgAAMoo/RFsaOsuxsSqo5sqg76p94vVeyQ7RtnE/FuO6MRkGMAAAkFbmor3ku2ofl5io7W5IunIs5KK9XwxgAAAgJS9Fe2utuoL+qvaxxETdiWjKQ7QLuWjvFwMYAAAYoXNHh1545lFVnr0sKXXR3lqr/cFjvqr28YmJVIdoF3rR3q+MA5gx5r9JulPSSWvtYve2aZKekdQs6YikDdbas8YYI+kJSZ+VdFHSfdbave7vuVfSX7t/7CPW2q1j+6UAAICxEFs5Vih90f5w8Lj6dNl31T7TIdrzPr5IG77xUD6+BePOyytg/13SP0p6Ou62r0t6yVr7mDHm6+7HfyXpDknz3X8+Iem7kj7hDmz/l6QVcgbbPcaY56y1Z8fqCwEAAKP37W99RRX7jkrKXLQPT5mqAzMbtb9hrjN8ZVG1T3eIdt+subp56SKtaf1qrr/8nMk4gFlr240xzQk33yVptfvzrZJeljOA3SXpaWutlfRrY0ydMabBve+L1tozkmSMeVHS7ZJ+POqvAAAAjFrnjg690LZRFZE+Sd6K9t31DUPjU5ZV+/SHaH++5FaOifxeAzbTWntckqy1x40xV7m3z5b0ftz9jrq3pbp9BGPM/ZLul6SmpiafDw8AAHgVO0Q7tnL0VbS3UQUU1TLbkbRqXxNel3LlWEyHaI+Vsb4I3yS5zaa5feSN1m6WtFmSVqxYkfQ+AABgbGzeuEm9e9tlNNT2yrZoH9CAbrPbtDK6U4s7DyWt2tdlqNqXUmLCC78D2AljTIP76leDpJPu7UclNcbdb46kHvf21Qm3v+zzcwMAgFGKr9rHhq/4laOUXdH+zt3bVf1aQKEj2VXti+0Q7bHidwB7TtK9kh5zf3w27vY/N8a0ybkIv9cd0v5V0qPGmKnu/T4jqTTf1gAAQIFLVrVPTExkW7Sv3h2SDShj1X7YyrFEqvZ+eMlQ/FjOq1f1xpijct7N+JikfzbGfFlSt6Q/cu/+L3ISFIflZChaJclae8YY818l/ca937diF+QDAIDcyVy1z/4Q7V/OCKlvldG+JpOxaq8SrNr7Yawt3MusVqxYYTs6OvL9MAAAKAleqvZDK8drPByiPfQOxzmRmboSrdGZ87doZpqqfeWMGfpCia4cjTF7rLUrvNyXEj4AACXOa9W+K3BUW5r6tWv+StnBqSxz0b4l/CntjXxWViZt1X761Y1qfaC1LFeOiRjAAAAoYV6r9kOJiSXO8OX5EO012hNZp1n9Yd2Y5HqvUjxIeywwgAEAUKK8rhxHc4h2T98CLTy3T6tP70p6vVcpHqQ9FhjAAAAoIZ07OnToN3v0zrGXNSF8XpK3qr2fQ7Rn9Z/Qqgsva9H5/QqkrNqX3kHaY4EBDACAEhEr2kvSBA1V7RfXXKuWamfl6KVqbxTVl+xmrYluU01bME3R/j3JHbzKuWrvBwMYAAAl4MlHHtSlzv2SnGEo4Fbtb5gYUl3Ie9V+eGIipK0tIQ9Fe/e3l2nV3g8GMAAAiphTtP+mQmecQ7QDwQYFqxZo3uQlaRIT6av2idd7ZS7aS8Ep0zW9Nqrb1t9SkomJscYABgBAkYov2sfWjfUT5quxKqDmyqDvqn3i9V4Zi/bNLbrnztlS80oSEx4xgAEAUITii/bDD9EOxN0ry5VjQmIiFL5dtSkO0aZoPzoMYAAAFJmReYmRh2hba9UV9L9yvNCdauU4VLTfUKJF+1xgAAMAoEh4KdqfML06HjijS7qiZxqNr5Xj0cg6ivbjjAEMAIAi4KVo3xU4qlcrDsrK+lo51oTXxiUmKNqPJwYwAAAKXDZF+/CUqTows1H7G+ZmXbWvOxFNeb0XRfuxxQAGAECB6tzRoefbHlNlpF+St6J9d33D0PiURdU+ceUYQ9F+fDCAAQBQgGJV+0olXzl6KdrLRhVUVPcNbBmxchxZtU+RmKBoPy4YwAAAKDDxVfusD9GWJGsV0IBus9u0MrpTy9veHbZyTF+1txTtc4ABDACAApFYtfd3iPbQtV537t6u6tcCCh4JZlG1N6qaMUN3k5gYVwxgAAAUgGRV+3mTWoYlJjKuHOOu9VrW9q72nQypt9ao3V05eqraz5mre/7yXhIT44wBDACAPNu8cZPO7W2XlKpqn+0h2kFtbanQ9pXO74+tHFNV7Vk55h4DGAAAeeKsHB9W6MwFSamr9n4P0W48e5XOfzRNx/pWs3IsMAxgAADkQfzKUUpdtfd7iHZL+FPaG/msrEyaqr2o2ucJAxgAADkWf5B2uqq930O0a8JrtCeyTrP6w7oxyfVeVO3zjwEMAIAcyqZqP2L48li07+lboIXn9mn16V1U7QsUAxgAADmQrmofWzlKmRITmYv2s/pPaNWFl7Xo/H4FqNoXLAYwAADGUeeODu39xS/1wXu/Tlm195KY8H6I9nuSO3hRtS9cDGAAAIyTJx95UJfe3i93C+i7ah8bvrI7RNv97SQmChIDGAAAYyyxaG8kBUZZtf9C244sD9GWzJSrNG1WnT5913pWjgWGAQwAgDE0omhftUDSJM2bfI2WTAq5q0F/Vfvn4or2A+HFHKJdxBjAAAAYI5s3blLv3nYZDRXtr6kKqSYoTQ6ZweHLb9Weon3pYAADAGCU4ov2seFr+Lox9tqU/6o9h2iXFgYwAABGwUvRXpJOZLFyTJaYSHaIdnzRnkO0iwsDGAAAPnkp2odNRIeDx3Ug2KN9vqr2a1Xb3ZB05UjRvngxgAEA4IPXov2rFQdlZX1X7etORJOsHCnaFzsGMAAAstC5o0MvPPOoKs9elpS5aH9g5vXqq6zSe9NnJQxfmav2qQ/Rpmhf7BjAAADwKLZyrJCPor3ko2qf/BDteR9fpA3feCjXXz7GEAMYAAAejOoQbcl51ctGFVRU92VYOSa73ksy6ps1VzcvXaQ1rV/N9ZePMcYABgBAGukO0c6maB/UFa22L2lldKfnlWPM0Mrx86wcSwQDGAAAKfziuz/Wvpd/NHiIdnDCcs2b1OK7aH/D7sOa8HpAbwWHqvY14XUpV45U7UsXAxgAAEl8b9PXdGHPO5KGqvYLJ4RUX+mMR36K9u1TKtR+izN4xar2dVTtyxIDGAAAcXoOdunpJ/5aoVOXJI1cOUqjK9rLiqo9GMAAAIgZcZB2bOUYl5g4YXr1VvCIr6L9rN7rdCZyk66cnTqiaj9s5UjVvuQxgAEAoOFV++Erx1hGIsuVY0JeojZ8mw5FPq2F595J+S5HqvblgwEMAFD2RiYmxm7l6OQl1upo34IU73J0hq/KGTO0gZVj2WAAAwCUrWyq9lua+rVr/krZwanMe9FeMlp4bl/Kqv30qxvV+kArK8cywgAGAChL2VftlzjDVxaHaA+EF2nVhXZNGrio5ovvxb3yxUHa5Y4BDABQVr7/nYd0/HeHNOmDfhn5qNpncYh2/LVezp8rcZA2JAYwAECZiC/aT3Zv81e1z/4Q7ZjYZfccpA0GMABAyUss2gcq5kiq1PV1n1RLtbNy9FK1N4rqS3az1kS3qaYtqH0nsyvaV0xt0NJ/c7NW3duah+8CCgkDGACgpD35yIO61Llf0lBeor4iqGlBaZabmMi+ah/S1paQtq8MULSHLwxgAICS1HOwS08//k2FzvRJSlg3SrEJKevExIRXK/SUu3KkaA+/GMAAACUnXdE+dq2XtVYnA+c8V+2Xtb07bOU4JbxWV8KLKdrDFwYwAEBJyVS0jw1f+4PHsjpIe2tLxbCVY22KlSNFe3jBAAYAKAk9B7v0wyceVvDUBUmpi/ZhE9Gh0HHtnNqXddU+/cqRoj28YwADABS92MoxllBNVrSPHaL9+tQB7Z/ZqAMNTYoa93d4TEwcTZKYoGgPPxjAAABFLX7lmKpoH/8Ox1fm36Coca/YyrJqnywxQdEefjCAAQCK1shDtLMs2vus2lO0x2gxgAEAik581V7yV7QPaEDLbIdqFdHK6E4fVXuK9vCPAQwAUFQSq/aJK0evRfvWgS26462XFDxvVPVawPPKcTAxcc31uudvv52PbwFKAAMYAKBoJFbts1o5SsOu9fpC2w5N3F2haEB66jPOq14zextU0f2ZlCtHqvYYKwxgAICCl7Zq7/kQ7eHXelX9KqgXlxq1X++EVRu7P6l3Lt6VduVI1R5jhQEMAFDQUlbtJ4Xc1WDmlWPSov2fOIOXsUZTwmt05szi9CtHqvYYQwxgAICCtXnjJp3b2y4psWrvjEfZH6I9vGhf07tIlyLLVcvKETnGAAYAKDjOyvFhhc6kr9pns3KMz0vMiszR73vv1Lm+uawckRcMYACAghK/cpTSV+29rByT5SUORdZKMmmq9qJqj3HFAAYAKBjZVu0zrhxTFO2X97+hioFLWn7uzbhXvqjaI3cYwAAABWHUVXt5Ldo/J6Po4LqRqj3ygQEMAJA3nTs69OLP/0kfXjqn6sglSclXjlKmxET6lWOyon38ulFU7ZFjoxrAjDH/u6Q/lfN/v52SWiU1SGqTNE3SXkn3WGsvG2OqJD0tabmk05LuttYeGc3nBwAUr1jRXpKqlXrl6KVq7/cQbet+Fqr2yDXfA5gxZrak/yRpobW2zxjzz5K+KOmzkh631rYZY74n6cuSvuv+eNZae40x5ouSNkq6e9RfAQCg6MQX7QPBBgWrFmje5CWjqtpnc4i2ldHlaTNVVzlBi5at0Kp7W3P+PUB5G+0KMiSp2hjzkaSJko5LWiPpj91f3yrpv8gZwO5yfy5JP5H0j8YYY621AgCUhWRF+6tq1qmxKqjmyqDvqr2fQ7TvZt2IPPI9gFlrjxlj/k5St6Q+Sf9T0h5JEWvtFfduRyXNdn8+W9L77u+9YozplTRd0im/jwEAUDwyFe2lUVTt/71Tta8Jr+MQbRSF0awgp8p5VetjkiKS/oekO5LcNfYKl0nza/F/7v2S7pekpqYmvw8PAFBANm/cpN697TJKXbS31qor6L9qHwrfrrrueor2KAqjWUGuk/R7a+0HkmSM+ZmkmyXVGWNC7qtgcyT1uPc/KqlR0lFjTEhSraQziX+otXazpM2StGLFCtaTAFDE4ov2seErWdHeWqv9wWO+qva1kUb19H5O1534kKI9isZoBrBuSZ80xkyUs4JcK6lD0g5JfyjnnZD3SnrWvf9z7se/cn99O9d/AUDp8lq0PxTsUZ8u+67aXwkv1vreNzhEG0VlNNeAvWaM+Ymc1MQVSb+V88rV/yepzRjziHvbD9zf8gNJPzTGHJbzytcXR/PAAQCFK5uifXjKVB2Y2aj9DXOd4SuLqn1td0PSlaMo2qPAmUJ+EWrFihW2o6Mj3w8DAJAFr0V7Z914jbrrG4bGp5Qrx+CwxMS5M+s080Q0ycqRoj3yxxizx1q7wst9KeEDAMZE544OvfDMo6o8e1lS+qL9lqZ+7Zq/UnZwKnMHLxtVQFEtsx0+q/YU7VEcGMAAAKMWWzlWyGvRfokzfMWtGwMa0G12m1ZGd2px56Gsq/YfTa7TvI8v0oZvPJSPbwGQFQYwAMCojOUh2nfu3q7q1wIKHfFetZeM+mbN1c1LF2lN61dz/vUDfjCAAQB86dzRoefbHlNlpF/SyMSEt6L98Hc4Vu8OyQaUddV+wz2fZ+WIosIABgDIWuwg7UoNVe0X11yrlmpn5eilaG8U1ZfsZq2JblNNW1AvXgmpb5XRviZD1R4ljwEMAJCV7236mi7seUfSUNX+hokh1YWGqvZeD9F2ivYhbW0JaftSZ0Cjao9ywAAGAPCk52CXnn7irxU6dUlS6qp9NodoxxftF4XnaW/wGl25OI+qPUoeAxgAIKOUB2knVO3fCh7xVbRvDn9Kr0fu0FX9J3UjVXuUAQYwAEBa8VX74QdpuwNWtivHJHmJzshaLTz3jlaf3kXVHmWBAQwAkNLIxMTYrRxjeYno2aladaFdi853KZCkal85Y4Y2sHJEiWEAAwAM07mjQy8++7gu9l3UhMhHknxW7TMcoh3LS6w+/ayMoklWjkbTr25U6wOtrBxRchjAAACD4teNE+S/ap9p5Zg6LyHe5YiywAAGAJA0fN0YCDYoWLVA8yYv8V21T7VyTFe0v1xXryu103THJ5o5SBsljQEMAMpcsqL9VTXr1FgVVHNl0HfV3ssh2jFWRpPmzNNNS2epecFyqvYoeQxgAFDGEov2I9eN/qr2+046K0eK9kByDGAAUKaefORBXercLyn1IdrWWnUFfVTtVwYo2gNpMIABQJnpOdilpx//pkJn+iSlPkTbWqv9wWOeExMTXq3QU+7KsTbSqJ7ez1G0B1JgAAOAMuKlaB82ER0OHlefLnuq2i9re3fYynFKeK2uhBdrPUV7ICUGMAAoE9kU7cNTpurAzEbtb5jrDF8pV45BbW2pGLZyrE2xcqRoDwxhAAOAEtdzsEs/fOJhBU9dkOStaN9d3zA0PnlITKRfOVK0BxIxgAFACYutHGPvaUxWtM94iLaNKqColtmOlImJo0kSE87YJYr2QBIMYABQouJXjqmK9pkO0Q5oQLfZbVoZ3anFnYeyqtobVo5ASgxgAFCCRh6iPTIx4bVof+fu7ap+LaDQkeyq9v2zGrX+1gUU7YEkGMAAoIQkq9qnOkTba9G+endINqCsqvaVsxq14Z4vULQHUmAAA4ASkblq77FoH/cOx1/OCKlvldG+JuPhIG2q9oBXDGAAUAK8VO0zXe8VG74S3+E4JzJTV6I1qui+JeXKkao9kB0GMAAoYl6r9ulXjolF+6F3ON548g+068z/IiuTduVI1R7IDgMYABSplFX7SSF3NZh55ZiuaN8c/pReidyhmf0ndCNVe2BMMYABQBHavHGTzu1tl5RYtXfGo2xWjqmK9p2RW7Xw3DtafXoXK0dgjDGAAUARcVaODyt0xnvVPtPKMVnR/qqzE7TqQrsWne9SgJUjMOYYwACgCHTu6FDnjhfVc2D34F/cvqr2CYmJVEX71af/VUbRpCtHqvbA6DGAAUCB2/oP39EHr26XuwX0XbUftnLMomjv/JuVIzCWGMAAoIDFF+2NpMAoq/aJK0cvRfvLdfW6UjtNd3yimao9MEYYwACgACUW7QPBBgWrFmje5CVaMrHCd9U+ceWYqWg/ac483bR0lpoXLKdqD4whBjAAKDCxQ7Tji/Ytk5pUEzSaHIqVvbKs2me1cqRoD4w3BjAAKCDf/tZXVLHvqKTkRXtpdFX7TCtHrvUCcoMBDAAKQOeODr3QtlEVkfRFe2utuoL+qvZeD9G+m0O0gXHHAAYAeRY7RLtC6Yv24cBZXdIVPdNosq7a14TXsXIECggDGADk0eaNm9S7t11G6Yv2r1YclJX1XbWv665n5QgUEAYwAMiD+KJ9bPhKV7Q/MPN69VVW6b3pszwnJmJV++tOfMgh2kCBYQADgByLP0RbyrJoL3lOTFwJL9Z6DtEGChIDGADkUCwxIfkv2stGFVRU9w1sSZmYqO1uSBlW/YiVI5B3DGAAkCPxVftkiQkvh2gHdUWr7UtaGd2p5W3vpkxMjFw5OsNX/6xGrb91AUV7IM8YwABgnHXu6NALzzyqyrOXJSVfOUpO0X5LU792zV8pOziVjVw3Ltl9WNWvBRQ8kjkxYd0/PZaY2EBiAigIDGAAMI5iK8dYYiLZynF40X6JM3yleYfjm8GQemuN2j1U7Q0rR6AgMYABwDjxunL0c4i2rFR/IfNB2qwcgcLEAAYAYyzxIO1UVXu/h2jP7L1OZyM36fd9CzJW7Vk5AoWJAQwAxlCsah87SDs4YbkW11yrlmpn5WhGrBxTHKJtN2tNdJtq2oLDivZTwmt1OLJOs/rDVO2BIsYABgBj5MlHHtSlzv2Shqr2N0wMqS40vGrv5RBt53qvkLa2hIYV7Y9FbtXCc+9o9eldVO2BIsYABgCj5FTtv6nQmeQHaUveEhOZivbRvmYtPLePqj1QAhjAAGAU4qv2ww7S9lK1z6Jov/LCQU0a6FLzxffihi+q9kCxYgADAB+ef+EZvdn+oia8N/Qux6GDtGPHBmW5ckwo2ofCt6s24RBt50+VqNoDxY0BDACyFMtLTHA/Ho+V44Xu1IdoO5/AqLJ+hjawcgSKEgMYAHiUWLQ3kgIZq/bZrxyPpi3aS8EpM3Xt3Cn67L3/jpUjUKQYwADAg5FF+xWaUTFdNZXT1VSVomqf5cqRoj1QPhjAACCDlEV7KTYhjapqH3+INkV7oDwwgAFACl7Or/sVAAAgAElEQVSK9tZ9dcpv1T7VIdoxFO2B0sQABgBJJBbtEw/RNjKy1upEoFdvBd/TL9NV7ZOsHA/NNqoJr0u5cqRoD5Q2BjAASJBYtE92iHbYRHQodFwHgj0ZExOJK8dY1b4uITFB0R4oHwxgAODKVLSPP0T7Z9M/0LG66bocXKA3m+anvd4rfuUYq9qnSkxQtAfKAwMYAMhb0T7+HY6vzL9FUTO0NEx2vdeytndHHKR9JbxY63vfSL1ypGgPlAUGMABlL5aYkHwU7a1Nfr3XtqC2tlQMO0g7sWpP0R4oXwxgAMpWz8Eu/fCJhxU8dUGS/6K97ICCsrovzUHayVeOzvBVOYOiPVBuGMAAlKXYyjGWUA0mKdpnOkTbKKrP2Wc10V7UQrvPR9XeaPrVjWp9oJWVI1BmGMAAlJ34lWOyxEQ2h2jfuW27Av1GFYeM3gpStQfgDQMYgLKSsmrvs2g/4dWQXlpi1H6Lc6E9VXsAXjCAASgL6ar2iYdoZ1O0fyqu7TUlvIaqPQBPGMAAlLxMVfsRh2hnUbQ/fHVAU3uv04eRmzRwto6qPQBPRjWAGWPqJD0labGcv2W+JOmApGckNUs6ImmDtfasMcZIekLSZyVdlHSftXbvaD4/AGTipWrv9XqvxHc4NpyfrYruteruW6CF596hag/As9G+AvaEpF9aa//QGFMpaaKkb0h6yVr7mDHm65K+LumvJN0hab77zyckfdf9EQDGXDZV+3SJiZRF+/BaHYyslWTSrhyp2gNIxvcAZoyZIulWSfdJkrX2sqTLxpi7JK1277ZV0styBrC7JD1trbWSfm2MqTPGNFhrj/t+9ACQRMqq/aSQuxrMvHL0UrRf3v+GKgYuafm5N+OGL6r2ADIbzStg8yR9IOn/NsYskbRH0l9Imhkbqqy1x40xV7n3ny3p/bjff9S9bdgAZoy5X9L9ktTU1DSKhwegHG3euEnn9rZLSqzaO+NRNivH9EX752QUHXzFi5UjgGyMZgALSVom6T9aa18zxjwhZ92Yiklymx1xg7WbJW2WpBUrVoz4dQBI1HOwSz/9/t/p7IdnVX32soz8V+3jV45eivZDUVXnM7ByBODFaAawo5KOWmtfcz/+iZwB7ERstWiMaZB0Mu7+jXG/f46knlF8fgAYXDdKzkWokr+qfbLERHzRPtkh2sOHL1aOALzzPYBZa8PGmPeNMddaaw9IWivpHfefeyU95v74rPtbnpP058aYNjkX3/dy/ReA0Ugs2gcq5mhaxXRdP2WRr6p9ssRETXitarsbkr7D0cjo0rSZqqucoEXLVmjVva25/hYAKFKjfRfkf5T0I/cdkL+T1CopIOmfjTFfltQt6Y/c+/6LnATFYTkZCv6mAuBbYtH+qpp1aqwKqrky6LtqH79yjC/apzpEm6I9AL9GNYBZa9+QtCLJL61Ncl8r6c9G8/kAoHNHh1545lFVnr0sKeFaL8Ve2PJXtY9fOSYr2scfok3RHsBoUMIHUDRiK8cKJS/aGxlZa3Ui4L9q7+UQ7XkfX6QN33goP98EACWBAQxAUfj2t76iin1HJaUu2ltrtT94zFfV3ssh2n2z5urmpYu0pvWruf7yAZQYBjAABa1zR4deaNuoikj6on3YRHQodFw7p/ZlXbVPtXKMGVo5fp6VI4AxwQAGoGDFDtGOrRzTFe1fnzqg/TMbdaChSVHjvgPSQ9W+Jrwu5cqRQ7QBjBcGMAAFafPGTerd2z4YVc1UtH9l/g2KGnd8yqJqX9ddzyHaAHKOAQxAQXEO0X5YoTMX/BftfVbtYzhEG8B4YwADUDDiD9GW/BXtAxrQMtuhWkW0MrrTc9WeQ7QB5BIDGICCkFi1T0xMeC3atw5s0R1vvaTgeaOq1wKeq/ZyExOsHAHkAgMYgLxLrNonS0xkc4j2xN0Vigakpz7jvWpfOWOGNrByBJAjDGAA8iZd1T62cpScov2Wpn7tmr9SdnAqS160r/pVUC8uNWq/3nmXo5eq/fSrG9X6QCsrRwA5wwAGIC8yVe0Vl5hwrvda4gxfmQ7R/hNn8DLWeKras3IEkA8MYAByzuvK0c8h2gtOXqe6U1N18cOWtFV7DtIGkE8MYABypnNHh55ve0yVkX5Jqav2fg/RnhtepY7IekVlPFTtOUgbQP4wgAHIiVjVvlJDVfvFNdeqpXroIO3hK8cUh2jbzVoT3aaatuCIov3bkbWa1R/WjekSE1TtARQABjAA4+7JRx7Upc79koaq9jdMDKkulLxqn+4QbadoH9LWltCwon00fK1WXWjXwvNdClK1B1DgGMAAjIueg1366ff/TpELEVVHLknyWbWXt6L96tPPyiga96qXg6o9gELEAAZgzMWK9pJUrYSDtLOo2qe73it90V6DV35RtQdQiBjAAIypxKJ9sGqB6qtvcA/SdgesbFeOscSEW7QPhW9XbYpDtGNFexsKafG181k5AihIDGAAxkxiXuKqmnVqrAqquTI4ZivHC92pDtF2hq8p9dN03dKr1XzjbbzLEUDBYgADMGrpivZGsRe2Rr9yPErRHkCJYAADMCqZivZGRtZa7Q/6XzlStAdQahjAAPjmpWjvDF/HfFXt4w/RpmgPoJQwgAHImteifdhEdCh0XDun9mVdtU91iHYMRXsAxYwBDEBWEov2yVaOsaL961MHtH9mow40NClq3EO246v2SVaOsap9qpUjRXsApYABDIBniUX7TIdovzL/BkWNOz4lud4rceUYq9rXpUhMULQHUCoYwABk1HOwS08//k2FzvRJ8nmIdpLrveJXjvFV+1QrR4r2AEoFAxiAtGJV+5BSF+0zHaId0ICW2Q7VKqKV0Z1a1vbusIO001ft3ZUjRXsAJYQBDEBKmzdu0rm97ZKGwqp+ivatA1t0x1svKXjeqOr1gLa2VAw7SDtV1Z6VI4BSxQAGYARn5fiwQmcuSBqbQ7Qn7q5QNCA99RlWjgDAAAZgmPiVozR8+PJ7iHbVr4J6calR+/VDK8f0VXtRtQdQ0hjAAAxKPEg7MTHh+xDtP3EGL2MNVXsAEAMYAJeXqn3K4StD0X7ByetUd2qqLn7YQtUeAMQABpS9dFX72MpRypCYSFO0nxtepY7IekVlqNoDgIsBDChjmar2XhITmQ7RfjuyVrP6w7oxXWKCqj2AMsMABpSpbKv26a73SnWIdvTsVK260K6F57sUJDEBAIMYwIAy0nOwS+/v69RL7T9XZc95ST6r9kpdtI8/RHv16WdlFI171csRWznezcoRQJliAAPKxK4f/1yvP/sDWavBleNg1X5SyB2SMq8c46/3Sla0HwgvTnGItgav/GLlCKDcMYABZWDzxk3q3ds+OBAFhlXtnaEom5XjHZ0vaeK2YNZFexsKafG181k5Aih7DGBACYsv2hu5r3pVLdC8yUuSXu/ldeUYf72Xl6J9bf00Xbf0ajXfeBsrRwAQAxhQspIV7a+qWafGqqCaK4O+q/aJ13slO0Q7vmjPIdoAMBIDGFCCMhftJd9V+7jERG13Q9KVI0V7AEiPAQwoMV6K9tZadQX9Ve1jiYm6E9EkK0eK9gDgBQMYUCI6d3TohWceVeXZy5JSF+2ttdofPOarah+fmEh+iDZFewDwggEMKAGxlWOF0hftDwePq0+XfVftMx2iPe/ji7ThGw/l41sAAEWFAQwoctkcoh2eMlUHZjZqf8NcZ/jKomqf7hDtvllzdfPSRVrT+tVcf/kAUJQYwIAile4Q7VRF++76hqHxKcuqffpDtD/PyhEAssAABhShxEO0fRXtbVQBRbXMdiSt2teE16VcOXKINgCMDgMYUGS+t+lrurDnHUlDba9si/YBDeg2u00rozu1uPNQ0qp9XYaqPYkJAPCPAQwoEj0Hu/T0E3+t0KlLkkauHKXsivZ37t6u6tcCCh3JrmpfNWOG7r7rE2r89Fdy+w0AgBLCAAYUgfiq/bCVY1xiItuiffXukGxAGav2w1aOVO0BYEwwgAEFLr5qP3zl6A5YPg7R/uWMkPpWGe1rMhmr9qJqDwBjjgEMKGAjExP+V46JeYk5kZm6Eq1RRfctaav2lTNmaAMrRwAYUwxgQAHyWrXvChzVlqZ+7Zq/UnZwKstctG8Jf0p7I5+VlUlbtZ9+daNaH2hl5QgAY4wBDCgwXqv2Q4mJJc7w5fkQ7TXaE1mnWf1h3Zjkei8O0gaA8ccABhSQbKr2fg/R7ulboIXn9mn16V1Jr/fiIG0AGH8MYEAB8FO193OI9qz+E1p14WUtOr9fgZRVew7SBoDxxgAG5FHnjg7t/cUv9cF7vx5WtV9cc61aqp2Vo5eqvVFUX7KbtSa6TTVtwTRF+/ckd/Ciag8A+cMABuTJk488qEtv75e7BRxMTNwwMaS6kPeq/fDEREhbW0Ieivbub6dqDwB5wQAG5FjPwS49/fg3FTrTJ8kZiAJjkJiY8GqFnvJctJeCU6Zrem1Ut62/hcQEAOQYAxiQQyOK9lULJE3SvMnX+K7aJx6i7alo39yie+6cLTWvJDEBAHnAAAbkSLKi/TVVIdUEpcmh2GX22VftEw/Rrk1xiDZFewAoHAxgwDjrOdilHz7xsIKnLkhK9g7H2Hjkv2qffuVI0R4ACg0DGDCOYivHWEI1VdH+RIZ3OWZKTByNrKNoDwBFhAEMGCfxK8d0RfvDweM6EOzRPi8rxxFV+7VxiQmK9gBQLBjAgHHgtWj/asVBWVnfVfu6E9GU13tRtAeAwsUABoyhbIv2B2Zer77KKr03fVbC8JW5ap+4coyhaA8AhY8BDBgjv/juj7Xv5R8NFu0TV45pi/bS8Kq9j5UjRXsAKB4MYMAYePKRB3Wpc78kH4doS86rXjaqoKK6L8nKcWZvgyq6P5Ny5UjRHgCKy6gHMGNMUFKHpGPW2juNMR+T1CZpmqS9ku6x1l42xlRJelrSckmnJd1trT0y2s8P5FNi1d7fIdpSUFe02r6kldGdI1aOjd2f1DsX70q7cqyaMUN3k5gAgKIxFq+A/YWkLklT3I83SnrcWttmjPmepC9L+q7741lr7TXGmC+697t7DD4/kBcjqvYTlmvepBYtmRRyV4OZD9GOv9brht2HNeH1gN4KDq/anzmzOP3Kcc5c3fOX95KYAIAiMqoBzBgzR9LnJP2tpAeMMUbSGkl/7N5lq6T/ImcAu8v9uST9RNI/GmOMte5LAEAR2bxxk87tbZc0VLVfOCGk+kq/h2gH1T6lQu23OINXpqo9K0cAKG6jfQXs7yV9TVKN+/F0SRFr7RX346OSZrs/ny3pfUmy1l4xxvS69z8V/wcaY+6XdL8kNTU1jfLhAWPLWTk+rNCZ5FV7aXRFe1l5OEiblSMAFDvfA5gx5k5JJ621e4wxq2M3J7mr9fBrQzdYu1nSZklasWIFr46hYMSvHKXkVftsDtFOvNZrVu91OhO5SUf7FqSp2ouqPQCUgNG8AnaLpM8bYz4raYKca8D+XlKdMSbkvgo2R1KPe/+jkholHTXGhCTVSjozis8P5IyXqn1WK8eEvERt+DYdinxas/rDVO0BoAz4HsCstQ9JekiS3FfAHrTW/okx5n9I+kM574S8V9Kz7m95zv34V+6vb+f6LxQDr1X7pMOXp6L92sFXvVaf3kXVHgDKwHh0wP5KUpsx5hFJv5X0A/f2H0j6oTHmsJxXvr44Dp8bGDPpqvbxB2mnv94rc9FeMlp4bh9VewAoI2MygFlrX5b0svvz30kacXGKtbZf0h+NxecDxlumqn3sIO10iQkvRftVF9o1aeCimi++Fzd8UbUHgFJHCR+Is+X7/1Un9r2tqhMfysh/1T42fHk9RNv5cyUSEwBQHhjAAA0v2k9wb/NbtY+/3iu7Q7QluSvHu1k5AkBJYwBD2Uss2gcq5kiq1MenfMJ31X5Z27vad3KoaF8TXpfxEO2KqQ1a+m9u1qp7W/PxbQAA5BADGMra5o2b1Lu3fXDdeFXNOtVXBDU5IDVVOQOWn6r91pYKbV8ZGCza11G0BwDEYQBDWYov2g9d6+WuG6XBbPBoqvYU7QEAqTCAoeykK9rHrvWy1upk4Jzvqv2U8FpdCS/W+t43OEQbADACAxjKSqaifWz42h885rtqXxNeq9ruhqQrR1G0BwCIAQxlxEvRPmwiOhQ6rp1T+7JeOcYnJkauHCnaAwCGMICh5HXu6NALzzyqyrOXJaUv2v9s+gfaP3OO9jfMzXrlmCwxMXSINkV7AMAQBjCUtNjKsULeivavzL9FUeNesZXFyjFVYiJ2iPa8jy/Shm88lIfvAACgEDGAoWSN+hDtLFaOqa736ps1VzcvXaQ1rV/N+dcPAChcDGAoOekO0fZatA9oQMvUoVob0croTs8rx5ihlePnWTkCAEZgAENJSTxEOzhhuRbXXKuW6qF3OXo5RLvVbtGa6DbJSjXPBLOu2nOINgAgHQYwlIzvbfqaLux5R9JQ1f6GiSHVhZzxKJui/Z3btivQH1TFYaOt80NU7QEAY4oBDEWv52CXnn7irxU6dUnSyJWjlH3RfsKrIb20xKj9ZudVL6r2AICxxACGopZ4kHZwwnLNm9QyLDFxwvRmXbR/yr3Wy1ijKeE1VO0BAGOKAQxFK75qH1s5LpwQUn2lO2Ble4h2Ql5iau91+jByk2rTvMuRqj0AwA8GMBSlkYmJ0a0c4/MSDednq6J7rbr7FqR4l6MzfFXOmKENrBwBAD4wgKGoZFO139LUr13zV8oOTmXeivYHI2slmbRV++lXN6r1gVZWjgAAXxjAUDSyrdrvmr/EGb6yLNov739DFQOXtPzcm3GvfA1V7Vk5AgBGiwEMRWHUVXt5Ldo/J6Po4Lox/novDtIGAIwVBjAUrM4dHXrx5/+kDy+dU3UkeWLCS9V+dIdoS+IgbQDAGGMAQ0GKFe0lqVrJV45eq/Zfspu1JrpNNW3ei/bGzUtIVO0BAGOPAQwF58lHHtSlzv2S3LZX1QLNm7zE+8pRGnm917aQtrZ4L9pfnjZTdZUTtGjZCq26tzXn3wMAQGljAEPB6DnYpacf/6ZCZ/okDbW9GquCaq4MZrFyTKzaVwyGVb0X7f+AvAQAYNwwgKEgpCvaG8Ve2Mq8coy/3mtZ27vDVo5Twmsp2gMACgIDGPIuU9HeyMhaq/3BbFaOQW1tqRi2cqxNsXKkaA8AyDUGMORNz8Eu/fCJhxU8dUFS6qK9M3wd81W1T79ypGgPAMgPBjDkRWzlGEuoJivanzC9OhTsUZ8uZ3WQdnxi4mjaxARFewBAfjCAIefiV46pivaxdziGp0zVgZmN2t8w1xm+sqzap0pMsHIEAOQTAxhyKtui/Svzb1DUuONT1lX75Nd7UbQHAOQbAxhyonNHh55ve0yVkX5J6Q/RTnmtl40qoKiW2Y6sqvYxlqI9AKBAMIBh3MWq9pXyeoj2yGu9AhrQbXabVkZ3anHnoaxWjoOJCYr2AIACwQCGcRVftc/6EG1p8Fqv1oEtunP3dlW/FlDoiPeVo5XRFa73AgAUGAYwjItkVfvsD9FOuNZrd0g2oKxWjk7VnsQEAKCwMIBhzKWs2k9yq/Y+ivYvXgmpb5XRviYzWLUfCC9Ov3Kkag8AKFAMYBhTmzdu0rm97ZISq/bOeJT9Idpu0X6pW8XPULVn5QgAKAYMYBgTzsrxYYXOpK/aZ7VyjMtLLAwv0G+DjbpycZ6Hg7RZOQIAChsDGEYtfuUopa7avxU84qto3xxerd9EPqOoTJqqvajaAwCKBgMYRiWbqr2nlWOSvERnZK1m9Yd1Y+8bVO0BACWBAQy+ZVu1HzZ8eSzaR89O1aoL7Vp4vktBqvYAgBLBAIasdO7o0Is//yd9eOmcqiOXJPms2mc4RDuWl1h9+lkZReNe9XJQtQcAFDMGMHgWv26slv+qvd9DtJ1/Ox9RtQcAFDMGMHjy7W99RRX7jkpy215VCzRv8hJfVfsvDWzO+hDtWF7ChkJafO18rvcCABQ1BjCk1bmjQy+0bVRFZKhof1XNOjVWBdVcGfRdtc/2EO1Jc+Zp8dJZal6wnJUjAKDoMYAhpdgh2hVKXrSX/FXt9510Vo6HZhvVhNdxiDYAoOwwgCGpzRs3qXdvu4xSF+2tteoK+qjarwwMFu3rKNoDAMoQAxiGiS/ax4avZEV7a632B4/5qtrXRhrV0/s5ivYAgLLFAIZBXov2h4I96tNlX1X7KeG1uhJerPVJoqocog0AKBcMYJCUXdE+PGWqDsxs1P6Guc7wlUXVvra7IenKURTtAQBlhAEMnov2sXVjd33D0PiURdW+7kQ0ycqRoj0AoPwwgJWxzh0deuGZR1V59rKk9EX7LU392jV/pezgVOYOXjaqgKJaZjsyVu2TH6JN0R4AUH4YwMpUbOUYS0xkLtovcYavuHVjQAO6zW7TyuhOLe48lFXVnkO0AQDljAGsDI3lIdp37t6u6tcCCh3xXrVn5QgAKHcMYGWkc0eHnm97TJWRfkkjExPeivbD3+FYvTskG1BWVXtWjgCAcscAViZiVftKDVXtF9dcq5ZqZ+XopWhvFNWX7GatiW5TTVtQL14JqW+V0b4mQ9UeAIAsMICVgScfeVCXOvdLGqra3zAxpLrQUNXe6yHaTtE+pK0tIW1f6gxoVO0BAMgOA1gJc6r231TozNBB2smq9tkcoh2fl1gUnqe9wWt05eI8qvYAAGSBAaxExVfthx2knVC1fyt4xFfRvjn8Kb0euUNX9Z/UjVTtAQDICgNYCYqv2g8/SNsdsLJdOSbkJULh29UZuVULz72j1ad3UbUHACBLDGAlZmRiYuxWjrWRRl3o/pyuOjtBqy60a9H5LgWSVO0rZ8zQBlaOAACkxABWAjp3dOjQb/bonePbNKFn5PVeyav22R+ifdTNS6w+/a8yiiZZORpNv7pRrQ+0snIEACANBrAit/UfvqMPXt0uWWmCPFbts1w5pivaO//mXY4AAGSDAayIxa8bjaTAKKv26Q7RTlW0v1xXryu103THJ5qp2gMA4BEDWBFKLNoHgg0KVi3QvMlLBocvP1V7L4dox1gZTZozTzctnaXmBcup2gMAkAUGsCKTWLSfMflTapnUpJqg0eSQcdeDHqv2SVaOFO0BABh/DGBFJLFon7hulLKr2ieuHCnaAwCQGwxgRSBT0T62brTWqivoPTERv3KsjTSqp/dzFO0BAMgB3wOYMaZR0tOSZkmKStpsrX3CGDNN0jOSmiUdkbTBWnvWGGMkPSHps5IuSrrPWrt3dA+/9Hkp2odNROHAWV3SFT3TaDImJpa1vat9J4dWjlPCa3UlvFjrKdoDAJATo3kF7Iqk/2yt3WuMqZG0xxjzoqT7JL1krX3MGPN1SV+X9FeS7pA03/3nE5K+6/6IFLwW7V+tOCgr6/Eg7aC2tlRo+8qhlWNtipUjRXsAAMaH7wHMWntc0nH35+eNMV2SZku6S9Jq925bJb0sZwC7S9LT1lor6dfGmDpjTIP75yBOz8Eu/fCJhxU8dUFS5qL9gZnXq6+ySu9Nn+U5MZF+5UjRHgCA8TQm14AZY5ol3SjpNUkzY0OVtfa4MeYq926zJb0f99uOurcNG8CMMfdLul+SmpqaxuLhFZXYyjGWUE1WtE95iLaUddU+fvhyxi5RtAcAYJyNegAzxkyW9FNJf2mtPWfi35KXcNckt9kRN1i7WdJmSVqxYsWIXy9l8SvHVEX7TO9wlI0qqKjuG9iSddXesHIEACAnRjWAGWMq5AxfP7LW/sy9+URstWiMaZB00r39qKTGuN8+R1LPaD5/KRl5iHb2Rfugrmi1fUkrozu1vO3drKv2/bMatf7WBRTtAQAYZ6N5F6SR9ANJXdba78T90nOS7pX0mPvjs3G3/7kxpk3Oxfe9XP81smqf7hBtL0X7JbsPq/q1gIJHsqvaV85q1IZ7vkDRHgCAHBjNK2C3SLpHUqcx5g33tm/IGbz+2RjzZUndkv7I/bV/kZOgOCwnQ9E6is9dEpJV7dMeop2paL8tqDeDIfXWGrV7Okibqj0AAPkwmndBvqLk13VJ0tok97eS/szv5ys1mar2oynay0r1F9KvHKnaAwCQP5Twc8xr1T79yjF10X7uuWt0/OzN+n3fgrQrR6r2AADkDwNYDqWs2k8KeT5EO13RviV8s/ZG7tTM/hPpV45U7QEAyCsGsBzZvHGTzu1tl5RYtXfGo2xWjqmK9nsit2rhuXe0+vQuVo4AABQwBrBx5qwcH1bojLeqvZeVY7KifbSvWQvP7WPlCABAEWAAG0fxK0cpy6q9x6L9lfBirbxwUJMGutR88b0kVXtWjgAAFBoGsHHw/AvP6M32FzXhPSes6rdqP2zlmKRoX9vdMOwdjs6fKlG1BwCgsDGAjbFH/+ZPVXEorAnuTOS3ap9q5RhftE+2bnQ+gVH/TKr2AAAUKgawMdK5o0MvPPOoKs9eluS8EhUYZdU+2coxWdE+ft1oplylqbPq9Om71lO1BwCgQDGAjYHYIdoViq0bV2hGxXTVVE5XU5XPqr2PQ7TnfXyRNnzjoXx8CwAAQBYYwEbp29/6iir2HZWUsG6UBs8JGE3V3ssh2n2z5urmpYu0pvWruf3iAQCALwxgPnXu6NALbRtVEUldtLfuq1N+q/beD9H+POtGAACKCAOYD7FDtGMrx2RFe2utTgR69VbwPf3SR9W+JryOQ7QBAChRDGBZ2rxxk3r3tssoddE+bCI6FDquA8Ee31X7uu56DtEGAKBEMYB5FF+0jw1fqYr2P5v+gY7VTdfl4AK92TTfc2IiVrW/7sSHFO0BAChhDGAeZFu0f2X+LYqaoaVhNlX79b1vcIg2AAAljgEsg1hiQvJRtLfWc2IisWofP3xRtAcAoLQwgKXx6N/8qSoPOscJpSvap3uHo+yAgrK6L+uqvTN8Vc6YoQ2sHAEAKCkMYEkkVu2TrRwlp2i/palfu+avlDXD1/79PGsAAAjUSURBVI1GUX3OPquJ9qIW2n0+qvZG069uVOsDrawcAQAoMQxgCUZW7UeuHIcX7Zc4w1eSdzjeuW27Av1GFYeM3gpmV7Vn5QgAQOliAIvjdeXo9RDtCa+G9NISo/ZbnLaXl6p9/ywO0QYAoNQxgMlZOT7f9pgqI/2SklftpewP0X7KXTcaazQlvMZj1f4LVO0BAChxZT+Axar2lRqq2i+uuVYt1c7K0YxYOaY4RNtu1proNtW0BYcV7af2LtSHkT/QwNk6qvYAAEBSmQ9gnTs69PbOHw2r2t8wMaS60FDV3ush2k7RPqStLaFhRfvuyK1aeO4dqvYAAGBQWQ9gz237vqps+qq910O0kxXto31z064cqdoDAFCeynoAO2X61ZShap/pEO1URfsb+4+rYuC4lp97M274omoPAADKfACrvdSkG6atVUt17NuQ5coxoWgfCt+u2u56rT79nIyig694UbUHAADxynoAm1+1XC1B59Wt0a4cL3SPPEQ7FlWlag8AAOKV9QC2pPp96dx0nQycG9XK8Wjaor0kqvYAACBOWQ9gkUu/U0/wKr3qc+WYqWh/adpM1VVO0KJlK7Tq3tb8faEAAKCglPUA1n5pri5POOipap/qEG2K9gAAIFtlPYAdnVylUzOmJhm+0q8cKdoDAIDRKOsB7P2GWu352IJhw1eqleOh2UY14XUpV44U7QEAgFdlPYAFq4/LmgUjhq/ElWOsal/XXU/RHgAAjFpZD2CLP/qN9upWfWStApLuc6/3il85xqr2iYmJGIr2AAAgW2U9gK169aDmNPwXvRNcpIV2n5a1vTvsIO1Y1X597xupV44U7QEAQJbKegA71lev5U/+Tjd+/F1VHDbaOr9i2EHatawcAQDAOCjrAey1RbV6biCsRd1W+242OjTbsHIEAADjrqwHsOs+mqRfzQnp0OwBGWs0JbwmY9Weoj0AABitsh7Allw8qtXv36K9Ez/Shx+2aOBsXcqqPYdoAwCAsVLWA9i/XLpRz9tbNXAxqIXn3qFqDwAAcqKsB7AD0yZp4IKhag8AAHKqrAewRfW/lz1erZWnd1O1BwAAOVPWA1jDmSpNP/2Ke5WXRGICAADkQlkPYKGTk3TJmsF3OMpdOd7NyhEAAIyjsh7AFt64VL8J/17WRiWxcgQAALlR1gPYqntbJUkH33xTLUuWDH4MAAAwnsp6AJOcIWxVvh8EAAAoK4F8PwAAAIBywwAGAACQYwxgAAAAOcYABgAAkGMMYAAAADnGAAYAAJBjDGAAAAA5xgAGAACQYwxgAAAAOcYABgAAkGMMYAAAADnGAAYAAJBjDGAAAAA5xgAGAACQYwxgAAAAOcYABgAAkGPGWpvvx5CSMeYDSe/l4FPVSzqVg8+D7PHcFDaen8LG81PYeH4Kl9/nZq61doaXOxb0AJYrxpgOa+2KfD8OjMRzU9h4fgobz09h4/kpXLl4blhBAgAA5BgDGAAAQI4xgDk25/sBICWem8LG81PYeH4KG89P4Rr354ZrwAAAAHKMV8AAAAByjAEMAAAgx8p6ADPG3G6MOWCMOWyM+Xq+H085Msb8N2PMSWPM23G3TTPGvGiMOeT+ONW93Rhj/sF9vt4yxizL3yMvfcaYRmPMDmNMlzFmnzHmL9zbeX4KgDFmgjHmdWPMm+7z87B7+8eMMa+5z88zxphK9/Yq9+PD7q835/PxlwtjTNAY81tjzAvuxzw/BcIYc8QY02mMecMY0+HelrO/38p2ADPGBCX9k6Q7JC2U9O+MMQvz+6jK0n+XdHvCbV+X9JK1dr6kl9yPJee5mu/+c7+k7+boMZarK5L+s7V2gaRPSvoz9/9HeH4KwyVJa6y1SyQtlXS7MeaTkjZKetx9fs5K+rJ7/y9LOmutvUbS4+79MP7+QlJX3Mc8P4XlNmvt0rjmV87+fivbAUzSTZIOW2t/Z629LKlN0l15fkxlx1rbLulMws13Sdrq/nyrpC/E3f60dfxaUp0xpiE3j7T8WGuPW2v3uj8/L+d/RGaL56cguN/nC+6HFe4/VtIaST9xb098fmLP208krTXGmBw93LJkjJkj6XOSnnI/NuL5KXQ5+/utnAew2ZLej/v4qHsb8m+mtfa45AwBkq5yb+c5yxN3HXKjpNfE81Mw3PXWG5JOSnpR0ruSItbaK+5d4p+DwefH/fVeSdNz+4jLzt9L+pqkqPvxdPH8FBIr6X8aY/YYY+53b8vZ32+h0fzmIpfsvyxochQ2nrM8MMZMlvRTSX9prT2X5j/KeX5yzFo7IGmpMaZO0s8lLUh2N/dHnp8cMsbcKemktXaPMWZ17OYkd+X5yZ9brLU9xpirJL1ojNmf5r5j/vyU8ytgRyU1xn08R1JPnh4LhjsRe2nX/fGkezvPWY4ZYyrkDF8/stb+zL2Z56fAWGsjkl6Wc61enTEm9h/X8c/B4PPj/nqtRq7/MXZukfR5Y8wROZe4rJHzihjPT4Gw1va4P56U8x8wNymHf7+V8wD2G0nz3XekVEr6oqTn8vyY4HhO0r3uz++V9Gzc7f/BfTfKJyX1xl4qxthzrz/5gaQua+134n6J56cAGGNmuK98yRhTLWmdnOv0dkj6Q/duic9P7Hn7Q0nbLSXucWOtfchaO8da2yznf1+22/+/nTtEqSCKwjj+/3gGQSxqNrgAV2AwGcQqCFrcg8kiCG8HZoug8IruwSoYNBtslhcF0zHcJ+7gzoP5/9JMO3Dg8t17z0zVKfZnKSRZS7L+9wwcAO90XN9G/Sf8JIe0HckEuK2q6cAljU6SB2Af2AK+gCvgCZgB28AncFxV80UguKF9NfkNnFfVyxB1j0GSPeAZeON/huWSNgdmfwaWZJc2JDyhbaZnVXWdZId24rIBvAJnVfWTZBW4o83yzYGTqvoYpvpxWVxBXlTVkf1ZDos+PC5eV4D7qpom2aTT+jbqACZJkjSEMV9BSpIkDcIAJkmS1JkBTJIkqTMDmCRJUmcGMEmSpM4MYJIkSZ0ZwCRJkjr7BVhu2SnY3c+5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab\n",
    "pylab.plot(covar, pheno,\".\")\n",
    "pylab.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
